Dynamical Thermalization of Disordered Nonlinear Lattices 
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We study numerically how the energy spreads over a finite disordered nonlinear one-dimensional 
lattice, where all linear modes are exponentially localized by disorder. We establish emergence 
of dynamical thermalization, characterized as an ergodic chaotic dynamical state with a Gibbs 
distribution over the modes. Our results show that the fraction of thermalizing modes is finite and 
grows with the nonlinearity strength. 
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The studies of ergodicity and dynamical thermalization 
in regular nonlinear lattices have a long history initiated 
by the Fermi-Pasta-Ulam problem m but they are still 
far from being complete (see, e.g., [4] for thermal trans- 
port in nonlinear chains and [3| for thermalization in 
a Bose-Hubbard model). In this letter, we study how 
the dynamical thermalization appears in nonlinear dis- 
ordered chains where all linear modes are exponentially 
localized. Such modes appear due to the Anderson local- 
ization, introduced in the context of electron transport in 
disordered solids 0, H, 0| and describing various physical 
situations like wave propagation in a random medium Q , 
localization of a Bosc-Einstcin condensate Q and quan- 
tum chaos 0]. 

Effects of nonlinearity on localization properties have 
attracted large interest recently. Indeed, nonlinearity 
naturally appears for localization of a Bose-Einstein con- 
densate, as its evolution is described by the nonlinear 
Gross-Pitaevskii equation An interplay of disorder, 
localization, and nonlinearity is also important in other 
physical systems like wave propagation in nonlinear dis- 
ordered media ll|, and chains of nonlinear oscillators 
with randomly distributed frequencies [l3l |. 

The main question here is whether the localization is 
destroyed by nonlinearity. It has been addressed recently 
using two physical setups. In refs. [3, [3it was demon- 
strated that an initially concentrated wavepacket spreads 
apparently indefinitely, although subdiffusively, in a dis- 
ordered nonlinear lattice. For a transmission through a 
nonlinear disordered layer , , chaotic destruction of 
localization leads to a drastically enhanced transparency. 

Here we study the thermalization properties of the dy- 
namics of a nonlinear disordered lattice - discrete An- 
derson nonlinear Schrodinger equation (DANSE). We de- 
scribe in details the features of the time-evolution of an 
initially localized excitation towards a statistical equi- 
librium in a finite lattice (we stress that this evolution 
is purely deterministic ~ and the relaxation to equilib- 
rium is due to determinsitc chaos.). Below we argue 
that a statistically stationary state is characterized by the 
Gibbs energy equipartition across the linear eigenmodes 



(Eq. (O) and call a relaxation to such an equilibrium 
state thermalization. Because thermalization is due to 
deterministic chaos, its rate heavily dependes on the sta- 
tistical properties of the chaos. As is typical for nonlinear 
Hamiltonian systems, depending on initial conditions one 
can obtain solutions belonging to a "chaotic sea" or to 
"regular islands". Moreover, one can expect the former 
to thermalize while the latter do not lead to thermaliza- 
tion. We numerically find non-thermalizing modes and 
characterize their dependence on the nonlinearity and the 
lattice length. We stress here that our analysis heav- 
ily relies on numerical simulations as analytic methods 
appear to be hardly applicable for disordered nonlinear 
systems. In numerics, a difference between thermalizing 
and non-thermalizing states (as well as between chaotic 
and non-chaotic states) is limited by the maximal inte- 
gration time: it might happen that the states which do 
not thermalize up to some time will thermalize in the 
future. There is no way to overcome this limitation in 
a simple way, because of a possibility for such slow pro- 
cesses like Arnold diffusion, characteristic time of which 
lies far beyond any computationally accessibility. Never- 
theless, peforming an analysis based on large but finite 
time scales, we can, on one hand, make predictions for 
experiments, and on the other hand, obtain a "coarse- 
grained" description of the dynamics. Accordingly, the 
results below should be understood as valid for available 
integration times, without a straightforward extrapola- 
tion for asymtotically large times. 

We describe a nonlinear disordered medium by the 
DANSE model: 
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Eni^n + /3| -0" I -011 + + i'n-1 , (1) 



where (3 characterizes nonlinearity, and the on-site en- 
ergies En (or frequencies) are independent random vari- 
ables distributed uniformly in the range —W/2 < En < 
W/2 (they are shifted in such a way that E ^ corre- 
sponds to the central energy of the band). We consider 
a finite lattice 1 < 71 < TV with periodic boundary condi- 
tions. Then DANSE is a classical dynamical system with 
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the Hamilton function 



n 

It describes recent experiments with nonlinear photonic 
lattices (cf. Eq. (1) in [l^l), where one follows, along a 
transversally disordered, finite nonlinear crystal, the evo- 
lution of a single-site or a single-mode initial state. This 
corresponds to the setup of our thermalization problem. 
Thus, the properties below can be observed experimen- 
tally as "thermalization of photons" , provided the crystal 
is long enough. In the context of many-particle quantum 
systems, Eq. ((2) is used as an effective mean-field Hamil- 
tonian of interacting bosons. 

For /3 = all eigenstates are exponentially localized 
with the localization length I w 96W^~^ (for weak dis- 
order) at the center of the energy band [6|. Below we 
mainly focus on the case of moderate disorder W ~ A, 
for which Z 6 at the center of the band and I Ki 2.5 
at £■ = ±2. For each particular realization of disorder a 
set of eigenergies £,„ and of corresponding eigenmodes 
^nm can be found. In this eigenmode representation 
Tpn = Cm'fnm the Hamiltoniau reads 



- P ^ VknjiCkCnCjC* , 
knji 



(3) 



with Em 1^™!^ = 1 ai^d Vjnm'mini{ ~ I arc the tran- 
sition matrix elements [3| • This representation is mostly 
suitable to characterize the spreading of the field over the 
lattice, since in this basis the transitions take place only 
due to nonlinearity. Also, the nonlinear correction to the 
energy is small (^ (3/1) for one excited mode. 

To study the dynamical thermalization in a lattice, we 
performed direct numerical simulation of DANSE ([T]) , us- 
ing mainly two methods: the unitary Crank-Nicholson 
operator splitting scheme with step At = 0.1 as described 
in [isl , and an 8-ordcr Rungc-Kutta integration with step 
At = 0.02; in both cases the total energy and the nor- 
malization have been preserved with high accuracy and 
both integration schemes gave similar results, for all lat- 
tice lengthes N used. Such a restriction of the accuracy 
check to the conserved quantities is suitable for chaotic 
systems. A comparison with other numerical methods 
for DANSE [l^ goes beyond the scope of this letter and 
will be performed in a longer publication. We started 
with two types of localized initial states: (A) one site 
seeded, i.e. |V'n(0)P = Snj and (B) one mode initially 
excited, i.e. |Cm(0)P ~ 6m,k- For different realizations 
of disorder, we seeded different possible sites/modes and 
followed the evolution of the field solving ([T]) up to times 
(in selected runs) ~ 10®. The level of spreading is char- 
acterized by the entropy of the mode distribution 



where overline means time averaging. For one excited 
mode 5* = while S* = IniV for a uniform distribution 
over all modes in a lattice of length TV. To give an im- 
pression of a time evolution of the thermalization process 
we show in Fig. [T] several representative time dependen- 
cies of the entropy One can see that for the setup 
(B) some modes remain localized during the complete in- 
tegration time (cf. [i^), while other after some transient 
time evolve to a state with large entropy. For setup (A) , 
the entropy grows in all cases with a tendency to satura- 
tion - some states seem to saturate at about S' « IniV, 
while others remain at values definitely smaller than In 
up to the maximal integration time. Especially the re- 
sults from (B) indicate a strong energy dependence of the 
spreading behavior, which is studied in this work. In our 
discussion below we focus therefore on the setup (B) as 
the mostly nontrivial one. 
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FIG. 1: (color online) Time evolution of entropy S ((4| in 
DANSE ig with iV = 32 and /3 = 1, for a particular re- 
alization of disorder and different initial states: bold black 
curves with markers - single-mode initial states (B) with en- 
ergies E = -0.34,0.76, -0.29,3.36, -0.5 (curves from top to 
bottom at t = 10*, two bottom cases are very close), solid 
red/gray curves - single-site initial states (A, ten randomly 
chosen states). The dashed line shows the level S = In 32. 
The time averaging has been performed over doubling time 
intervals (between successive markers). 



To derive an approximate expression for the statis- 
tically stationary distribution pm, we mention that it 
should satisfy = 1 and E = J^Pm^m^ where, in 

view of discussion above, we have neglected the nonlin- 
ear contribution to the energy. Then the condition of 
maximal entropy (U) leads, after a standard calculation, 
to a Gibbs distribution: 



Pn 



= z- 



exp(-e„i/T), Z = ^exp(-e,„/r). (5) 



Here T is an effective "temperature" of the system: it 
has no meaning as a physical temperature, but serves as a 
parameter characterizing the distribution; it is a function 
of the total energy E of the state and of the realization 
of disorder. The entropy and the energy are related to 
each other via usual expressions, e.g. |2l| : 



5* = - Pm In Pn 



Pn 



\c„ 



(4) 



E = T^dliiZ/dT , S ^ E/T + \nZ 



(6) 
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This value of entropy yields the maximal possible 
equipartition for the given initial energy, and the values 
of Fig. [T] obtained via a numerical simulation of the dis- 
ordered nonlinear lattice should be compared with these 
values from the Gibbs distribution. Because we have 
anyhow neglected the effects of nonlinearity in the the- 
oretical value of the entropy, we adopt other simplifi- 
cations: approximate the density of states of the disor- 
dered system as a constant in an interval — A < e < A 
and consider the energy eigenvalues in a particular 
realization of disorder as independent random variables 
distributed according to this density. Furthermore, we 
assume the variations of the partition sum to be small 
and use (InZ) w ln(Z), where brackets denote averaging 
over disorder realizations. In this simplest approximation 
we obtain 

(InZ) wlniV + lnsinh(A/T) -ln(A/T). (7) 

At = 4 we have A « 3 (see Figs. [313] below) and 
this theory gives the dependence S{E) within a few 
percent accuracy compared to S averaged over disorder 
within Gibbs computations with exact numerical values 
e„i- This justifies, for the parameters used, the approxi- 
mation above. We note that T = +0, ±00, —0 correspond 
to E = —A, 0, -l-A respectively (as in the standard two- 
level problem, see related discussion in [2l[). 
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FIG. 2: (color online) Left: time and disorder averaged prob- 
ability {pmijn')) in mode m for initial state in mode m' . 
Right: theoretical values according to the Gibbs distribution 
dSl). Here N = 2,2, (5 ^ I, Ni = 15. 

We compare in Fig. [2] the Gibbs distribution ([5]) with 
the results of direct numerical simulations of DANSE us- 
ing Nd disorder realizations. Here we present the values 
{pm) averaged over time and over different realization 
of disorder, in dependence of the number of the initially 
seeded mode m! . The modes have been ordered according 
to their energy, so that to = 1 corresponds to the maxi- 
mal energy. One can see a good correspondence between 
the numerics and the simple theory ([5]) in the sense that 
states at the band edges remain localized, while states 
in the center generally spread. However, there is one 
crucial discrepancy: the peaks on the diagonal m = to' 
indicate that there are cases when there is no thermaliza- 
tion within our simulation time and the energy remains 
in the initially seeded mode. 
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FIG. 3: (color online) Left panel: Final entropies Q after 
an evolution during time interval 10^, averaged over a time 
interval of 10®. The states evolving from initial modes in the 
middle of the band (see text) are marked with black circles, 
while those at the edges of the band are marked by the red 
(gray) pluses. The curve shows the approximate theory (0. 
Right panel: Lyapunov exponents " (averaged over a time 
interval 10®) vs. entropy for the same sets. Here A'^ — 32, 
/3 = l,iVd = 7. 



To characterize thermalized and non-thermalized cases 
quantitatively, we compare in Fig. [3] numerical values for 
S{E) according to Eq. (jl]) with the theoretical Gibbs 
computation given by Eas. (|5l6l7p . Clearly, the Gibbs 
theory gives a satisfactory global description of numeri- 
cal data. The nonthermalized modes in this presentation 
are those at the bottom of the graph; these states are 
absent for the setup (A) where initial sites are seeded. 
(Again, as discussed above, "nonthermalized" should be 
understood as "nonthcrmilized whithin the integration 
time"). 

It appears appropriate to discuss the dynamics of the 
modes in the middle of the energy band (|em| < 2) and 
at the edges (|em| > 2) separately. For the modes in 
the middle of the band, the maximal entropy according 
to ^ is close to InA^, and one clearly distinguishes the 
thermalized modes and those that remain localized, as 
those reaching the maximal entropy and those remain- 
ing at the level < 1, correspondingly. Thermalization 
is associated with the chaotic dynamics of the DANSE 
lattice. To demonstrate this, we calculated the largest 
Lyapunov exponents A shown in Fig. [3] (right panel). All 
modes with S < \, i.e. those that do not thermalize, 
have nearly vanishing A, while for the thermalized states 
[S > 2) the positive values of A clearly indicate chaos. 

The above distinction between thermalized and non- 
thermalized states is less evident for modes at the band 
edges (shown by red pluses in Fig. [3]) . Here already the 
theoretical value of entropy given by Eqs. (|5|6|7p is rather 
small. Hence, the spreading can go over a few "available" 
modes only. Nevertheless, also here one can see from 
Fig. [3] a clear correlation between the entropy and the 
Lyapunov exponent. Moreover, in several cases the Lya- 
punov exponent at the edge of the spectrum is definitely 
larger than in the middle. This happens because the en- 
ergy spreads over a small number of modes, hence the 
effective nonlinearity is larger due to larger amplitueds 
of each mode, and therefore chaos is stronger. 



Above, we did not account for a spatial organization of 
the mode structure. The latter is less important for the 
modes in the middle of the band, where one can always 
expect to find neighbors with a close energy. Contrary 
to this, for the energies at the edges the issue of spatial 
distance becomes essential. Indeed, since here the ther- 
malization is possible only over a few modes, it is impor- 
tant whether these modes are spatially separated or not. 
For linear eigenmodes m and m' the natural measure of 
this separation is the coupling matrix element Vm'm'm'm 
according to ([3]). It is exponentially small for spatially 
separated modes due of their localization. One can ex- 
pect that a mode at the edge of the spectrum will be 
thermalized only if the coupling V to other few modes in 
the lattice with a close energy is large, what is a rather 
rare event. 
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FIG. 4: (color online) Dependence of entropy S on energy 
E as in Fig. [2] but for N = 64, Nd = 18, and two values of 
nonlinearity: (a) /3 = 0.5, (b) P = 2. Averaging have been 
performed over the time interval 10^ after an initial evolution 
during time 10^; for small (3 still longer times are needed to 
reach thermalized state with maximal S at given E. 

Finally, we discuss how the thermalization properties 
depend on the nonlinearity constant (3. In Fig. [4] we show 
the dependence S{E) for different nonlinearities /3. For 
(3 = 0.5 a large portion of the initial states remains non- 
thermalized, while for = 2 all states are thermalized 
(at least in the sense that their entropy is close to the 
maximal possible one, as discussed above this is a good 
criterion in the middle of the band). To determine how 
the fraction of thermalized states depends on nonlinearity 
/3 we use the following procedure. For the initial modes 
in the middle of the band (i.e. for \E\ < 2) we classified 
those that reach more than the half of the maximal en- 
tropy (i.e. the level ln(iV)/2) as thermalized, and those 
that remain below this level as non-thermalized. The 
fraction fth of the thermalized modes, shown in Fig. [5l 
monotonously increases with At fixed /? the numerical 
data indicate saturation of fth at large TV, but more de- 
tailed checks at larger sizes and longer times are required. 
For example, recent results on self-induced transparency 
of a disordered nonlinear layer [l3| show decrease of crit- 
ical f3 with lattice size for iV < 32. 

The properties of thermalization described above can 
be incorporated in a general framework of nonlinear dy- 
namics as follows. One can expect, based on general 
KAM arguments, that for small nonlinearity regular. 
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FIG. 5: (color online) (a) Fraction of thermalized (after time 
10®) modes fth from the middle of the band as a function of 
nonlinearity (3 for A'^ = 16 (circles), 32 (bold line), and 64 
(pluses). Diamonds show data for t = 10^ and N = 32. The 
dotted line shows the fraction of the bifurcated breathers fb 
according to panel (b). Panel (b): the bifurcation values 13c 
for different modes vs their linear energies for A = 32. To 
all modes with /3c > 2 we have attributed Pc = 2, this set 
looks like two vertical "lines" at /3 = 2 on panel (b). 



non-ergodic regimes predominate, while for large values 
of (3 stable solutions are destroyed and a chaotic ergodic 
state establishes in the lattice. While it is hard to char- 
acterize this transition via a general analysis of the dy- 
namics in a high-dimensional phase space, it is possible 
to follow the evolution, as nonlinearity increases, of spe- 
cial resonant modes that stem from linear ones. Looking 
for solutions of ([1]) in the form 'ipnit) = (^ne"**^*, we arrive 
at a nonlinear eigenvalue problem e(/)„ = + I3(t>^ + 

<j)n^i + 4>n+i which, of course, at /? = yields linear 
frequencies and modes. Starting from these modes, we 
followed these solutions to larger nonlinearities using a 
numerical continuation, and in this way obtained nonlin- 
ear resonant modes - "breathers" (cf. 13, 23|). Worth- 
noting, these modes change in the regions where the field 
is large, while in the tails they follow linear solutions, in 
accordance with psj . Moreover, we performed numerical 
stability analysis of these breathers and found that they 
bifurcate at some critical value of nonlinearity Pc- The 
values of (3c for an ensemble of realizations of random po- 
tentials are shown in Fig. O). Additionally, we show in 
Fig. [5^ a cumulative distribution of /3c for the same range 
of eigencnergies |e„| < 2 that is used for the other curves 
plotted. First of all, note the similar global behavior of 
/th and /b which makes us believe that the bifurcations 
of stable resonant modes is indeed the mechanism of the 
/3-dependence of thermalization. However, the curves do 
not coincide because (3c is defined as the value of the first 
bifurcation, which may not immediately lead to chaos 
but may be the first one in a series of transitions to more 
irregularity. Strictly speaking, /b should be an upper 
bound for /th, which is seen in Fig. [5^. The increase of 
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/th from t = 10^ to 10^ shows that it hasn't saturated 
yet, but the saturation curve must he below fh- Re- 
markably, we have found that the breathers at the edges 
of the band, i.e. for |e„| > 2, are extremely stable: most 
of them remain stable up to large values of /3 « 5. This 
corresponds to the numerical observation of the strong 
suppression of the thermalization for these modes. We 
emphasize here that because of the nonlinearity of the 
system the superposition principle does not hold. This 
means that to observe a stable breather mode one has to 
prepare initial conditions mostly close to this solution - 
what is achieved here by choosing the initial conditions 
as a pure linear eigcnmode (case B above). When one 
initially seeds one site, as in case A (or uses other ini- 
tial conditions not close to a breather), then this initial 
condition does not produce a breather because the lat- 
ter typically does not survive nonlinear interaction with 
other components of the solution. If, for example, one 
starts with an excitation of two modes which are both 
stable at some value of j3, one might still see fast ther- 
malization, because a superposition of two breathers is 
not a breather. 

Our main conclusion is that the maximally thcrmalizcd 
state in a disordered nonlinear lattice ([T|), that emerges as 
a result of chaotic dynamics, is described by the Gibbs 
distribution over the linear modes, with some effective 
temperature depending on the initial excitation. Not all 
modes lead to thermalization, some fraction of them re- 
mains localized, but this fraction decreases with nonlin- 
earity. We found that this can be explained by the disap- 
pearance (via bifurcations) as the nonlinearity increases, 
of stable resonant modes - breathers - stemming from 
linear eigenstates. Further studies arc still required to 
establish the properties of this thermalization in depen- 
dence on the nonlinearity strength, disorder and lattice 
size. 
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